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Abstract 

A fast multigrid solver for the steady incompress- 
ible Navier Stokes equations is presented. Unlike 
time-marching schemes, this approach uses relaxation 
of the steady equations. Application of this method 
results in a discretization that correctly distinguishes 
between the advection and elliptic parts of the oper- 
ator, allowing efficient smoothers to be constructed. 
Numerical solutions are shown for Bow over a Bat 
plate and a Karman-Trefftz airfoil. Using collective 
Gauss-Seidel line relaxation in both the vertical and 
horizontal directions, multigrid convergence behavior 
approaching that of O(N) methods is achieved. The 
computational efficiency of the numerical scheme is 
compared with that of a Runge-Kutta based multi- 
grid method. 

Introduction 

One of the critical needs in computational fluid dynam- 
ics is faster flow solvers. Multigrid is a well known method 
of convergence acceleration that is widely used in Euler 
and Reynolds-averaged Navier-Stokes codes. These appli- 
cations of multigrid generally are based on the unsteady 
equations using some temporal integrator as the smoother, 
combined with a full-approximation scheme (FAS) multi- 
grid iteration. A common approach is one originally pro- 
posed by Jameson. 1 * Starting with the unsteady equations, 
a finite-volume spatial discretization with explicit artificial 
viscosity is combined with a Runge-Kutta (R-K) time in- 
tegration as a smoother. An alternative approach 2-4 is to 
use upwind- differencing and implicit time integration as 
the smoother. However, these approaches have resulted in 
poor multigrid efficiency. When applied to high Reynolds 
number flows over complex geometries, convergence rates 
are often worse than 0.99. There is clearly a need to de- 
velop substantially more efficient multigrid solvers. 

According to Brandt, 5 one of the major obstacles to 
achieving better multigrid performance for advection dom- 
inated flows is that the coarse grid provides only a fraction 
of the needed correction for smooth error components. 
This obstacle can be removed by designing a solver that ef- 
fectively distinguishes between the elliptic, parabolic, and 
hyperbolic (advection) factors of the system and treats 
each one appropriately. For instance, advection can be 
treated by space marching, while elliptic factors can be 
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treated by multigrid. The efficiency of such an algorithm 
will be essentially identical to that of the solver for the el- 
liptic factor only, and thereby attain so-called “textbook” 
multigrid efficiency. 

Relatively little research has been done in the area of 
multigrid algorithms for the Navier-Stokes equations based 
upon factorizable discrete schemes. Brandt 5 has presented 
an approach called “distributive relaxation” by which one 
can construct smoothers that effectively distinguish be- 
tween the different factors of the operator. Using this 
approach, Brandt and Yavneh have demonstrated textbook 
multigrid efficiency for the incompressible Navier-Stokes 
equations. 6 Their results are for a simple geometry and a 
Cartesian grid, using a staggered-grid discretization of the 
equations. 

Recently, Thomas, Diskin, and Brandt ' achieved text- 
book mult igrid efficiency for high Reynolds number incom- 
pressible wake and boundary layer flows associated with a 
flat plate. Their scheme uses a staggered grid approach 
with distributed relaxation and defect correction. With the 
distributive relaxation the system of equations was decom- 
posed (i.e., factorized) everywhere, except near boundaries 
where they remained coupled. In all calculations Cartesian 
grids were used. 

Sidilkover and Asher 11 introduced a fast multigrid solver 
for the incompressible Navier-Stokes equations that does 
not require a staggered grid arrangement of flow variables. 
With a pressure Poisson formulation as a foundation, a fac- 
torizable discrete scheme is derived. In this work also only 
simple model problems were solved using uniform Carte- 
sian grids. 

In this paper, an alternative to distributive relaxation 
is presented. It is a generalization of the approach of 
Sidilkover and Ascher 11 and a continuation of the work 
of Roberts and Swanson. 10 This approach can be classi- 
fied as a method of the Weighted Gauss-Seidel type. 5 A 
conventional vertex-based finite-volume or finite-difference 
discretization of the primitive variables is used, avoiding 
the need for staggered grids. This simplifies the restriction 
and prolongation operations, because the same operator 
can be used for all variables. A projection operator is ap- 
plied to the system of equations, resulting in a Poisson 
equation for the pressure. In the case of the incompressible 
Navier-Stokes equations the physical diffusion terms in the 
pressure equation require approximation only at bound- 
aries. A suitable boundary condition for general geometries 
can be derived. The Poisson equation for the pressure may 
be treated by Gauss-Seidel relaxation, while the advection 
terms of the momentum equation are treated by space- 
marching. Because the elliptic and advection parts of the 
system are decoupled, ideal multigrid efficiency is possible. 
With the present scheme line relaxation is used to treat 


1 

American Institute of Aeronautics and Astronautics Paper 



AIAA-2001-2574 


the anisotropy of mesh aspect ratio that generally occurs 
in resolving viscous flows. 

Mathematical Formulation 

The incompressible Navier-Stokes equations in primitive 
variables are 

UU X + VUy + P X = 1>{U XX + Uyy), 

UV X + VVy + Py = l>(v xx + Vyy), 

n x +Vy = 0, 

where u and v are the components of the velocity in 
the x and y directions, respectively, p is the pressure, and 
the subscripts denote partial differentiation. The density 
is taken to be one. The advection-diffusion operator is de- 
fined by 

Qv = Q ~ vA, (1) 

where Q = ud x + vd y , A is the Laplacian operator, 
and d x , d y are the partial differentiation operators. The 
coefficient v is the reciprocal of the Reynolds number, 
v = 1 /Re. The incompressible Navier-Stokes equations 
may be written as 

/ Qv 0 d x \ /u\ 

Lq = I 0 Q v <) v ] ( v 1=0. (2) 

\d x d v 0 J \pj 

Introducing the adjoint to Q, defined by 

Qt(f) = -dx(uf) - dy(vf) + A(vf), (3) 

a projection operator P is defined: 



Fig. 1 Primary and dual cells on quadrilateral 
grid. 

Discretization 

The first step in approximating L is to discretize the 
Navier-Stokes equations (2). Consider a typical grid vertex 
as shown in Fig. 1. The momentum equations are dis- 
cretized using a second-order accurate upwind difference 
stencil for the advection part of operator (1), and central 
differencing for the physical diffusion contribution. Central 
differencing is also used for the pressure gradient term. 

Consider the discretization of the term ud x u in the x- 
momentum equation, assuming a Cartesian grid for sim- 
plicity. The second-order upwind-difference discretization 
of this term may be written as 


// ° °\ 

P = 0 I 0 . (4) 

\d* 9, Ql) 

Applying the projection operator to the Navier-Stokes 
equations yields 


ud x u | 


so 



/2 ... ("i... 



u i- l,j) 

/2,j(n*— i Ui— 2 j ) (baj 


when u > 0, and 


/Qv 0 d x \ / u\ 

Lq = PLq = I 0 Q v I v + s.p.t., (5) 

\ 0 0 A) \p) 

The matrix operator on the right-hand side consists of the 
principal part of L, and “s.p.t.” are the subprincipal terms, 
in the terminology of Brandt.® These terms arise because 
the coefficients u and v in the operators Q„ and Ql are 
not constant. It is important to note that the subprincipal 
terms can be ignored for the purpose of constructing a 
relaxation scheme. 

The operator on the left-hand side of Eq. (5) is upper 
triangular. The pressure satisfies a Poisson equation for 
which a conventional relaxation method, such as Gauss- 
Seidel, can be applied. Upwind differencing of the advec- 
tion operator in the momentum equations allows down- 
stream relaxation to be used. The strategy used to relax 
the system is to first update the pressure. The pressure 
update contributes to the velocity update through the gra- 
dient terms in the right-hand column of the operator in 
Eq. (5). Finally, the velocity components are updated by 
relaxing in the streamwise direction. 


'U() x v | sd — 9/j 


“ ^■« I + 3/ 2j («, + 2, 


«.+ ij) (6b) 


when u < 0. Here, the superscript h denotes the discrete 
approximation to the corresponding differential operator, 
and the average velocity components Ui±!/2,j are 


t, s+l/2,j — 2 ( U '+ "b J')i 

1 / ^ 
u i — 1/2, j — + u i— l,j)j 


Analogous expressions may be written for the vdy opera- 
tor. 

The projection operator P is applied to the discrete 
equations to obtain the residual for the pressure Poisson 
equation of Eq. (5). Letting R Pij be the pressure equation 
residual at vertex (i. j), the application of P can be written 
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in integral form, 



where A tJ is the area of the control volume centered on 
vertex ( i,j ). This control volume is the area bounded by 
the dual grid cell, shown as the dashed lines in Fig. 1. 
The discretization of the pressure equation is done by first 
discretizing the boxed terms of Eq. (7) on the edges of the 
dual grid cell. The evaluation points are shown as solid 
squares. To evaluate the gradient of pressure p at the dual 
grid face (i — 1/2, j) the partial derivatives are written as 

Px = Piix +Pr,’lx, Py = P(£y + PyPy, (8) 

where { and p are the generalized coordinates correspond- 
ing to the i and j directions on the grid respectively. The 
derivatives d^p and d v p are approximated at the face center 
by 

diP\ = Pi,] ~ Pi-i, j, 

dyP>\ = t (Phj+i +p<- 1.J + 1 — Pi, j-1 — Pi- i,j — i ) , 

h-l/2 ,1 4 

0) 

with similar expressions for the faces (/ | 1/2../}; ( i , j-f-1/2), 
and (i, j — 1/2). The gradients of u and v are found the 
same way. The grid metric terms £x, £y< Px, and p y are also 
evaluated on the dual grid face centers. These expressions 
are used in the boxed terms in Eq. (7). The integral (7) 
is then evaluated to get the pressure equation residual at 
the vertex (i, j), taking the boxed terms to be constant over 
each face of the dual grid cell. For a uniform Cart esian grid, 
the principal part of the resulting stencil is a conventional 
five-point approximation to the Laplacian operating on the 
pressure . 

With the velocity gradients known at the faces of the 
dual grid cell of Fig. 1, one can easily determine the physi- 
cal viscous terms in the momentum equations. By applying 
Green’s theorem, as in Eq. (7), to the dual grid cell, the 
Laplacian in A u and Av is approximated by 



8A ii 


Viscous terms also appear in the pressure equation, 
which can be written as 

A p = -dx[uu x + vuy)] - d y [uv x + VVy)] 

+ i^&k/Au) + c^(Au)], (11) 

where the terms on the right side of (11) are subprin- 
cipal terms with respect to relaxation. If the functions 
u(x,y) and v(x,y), as well as their spatial derivatives, are 
continuous, then the order of differentiation does not mat- 
ter. Consider a two-dimensional flow domain where these 


functions and their derivatives are continuous. Then, in 
the interior of the domain, the derivatives of the viscous 
contributions cancel with the application of the continuity 
equation. 

At a solid boundary the viscous terms do not cancel. 
The inviscid subprincipal terms of Eq. (11) vanish when 
the no-slip conditions are applied, as in viscous flows. The 
pressure equation reduces to 

A p = v[d x ( Au) + dy(Av)], (12) 

In the discrete problem ( 12 ) is approximated for a half dual 
cell at the boundary. From the momentum equations, 

Px\,,i = u(Au)i,i, ^ 

p v \i,\ = v(AD)j,l, 

where j = 1 corresponds to the surface boundary. The 
discrete form of Eq. (12) provides the necessary boundary 
condition for pressure. The present approach for treating 
the viscous terms in the pressure equation is similar to that 
of Sidilkover and Asher. 11 

Solution Procedure 

Rather than discretizing Eq. (5) directly, the Erst step 
of numerical solution procedure is to discretize Eq. (2). 
The relaxation scheme is constructed by applying the pro- 
jection operator P at the discrete level rather than the 
differential level. A sequence of grids Gk , ( > K i , . . . , Go is 
used, where Gk is the finest and Go the coarsest. Let L* 
be the discrete approximation to the operator L and q* 
be the solution on the fc-th grid. This system has the 
form L;. ([;, = f*, where the entries of L* are 3x3 block 
matrices which operate on the unknowns («, v,p) T at each 
grid vertex. A general iteration scheme is constructed by 
writing the operator L k as L/, = M* — NT, where the 
splitting is chosen such that Mfc is easily inverted. Lexico- 
graphic Gauss-Seidel is obtained by taking Mfc to be the 
block lower-triangular matrix resulting from ignoring the 
blocks above the diagonal of L*. A further simplification 
is obtained if the diagonal blocks of Mfc contain only those 
entries corresponding to the principal part of the operator. 
Because the operator in Eq. (5) is upper triangular, the di- 
agonal blocks will then be 3 x 3 upper triangular matrices. 

Letting q/ be the n-th iterate of the solution, the itera- 
tion is 

M; «i; " = f; + N. q}' . 

Because the operator L;, : is nonlinear, Mt and NT will be 
functions of q£ and q/ +1 . Letting i5q/ = q/ +1 — q”, the 
iteration may be rewritten as 

M*iq2=f*-L*qg. (14) 

Because Mfc is block lower-triangular, <S([" is found by for- 
ward substitution and inverting a 3 x 3 diagonal block at 
each vertex. The diagonal blocks are upper triangular and 
are easily inverted. 

If upwind differences are used for the advection opera- 
tor (1) and the grid points are ordered in the flow direction, 
then the 3x3 blocks of N* will have zeroes in the first two 
rows. In this case, lexicographic Gauss-Seidel relaxation is 
equivalent to space-marching of the advection terms. The 
advected error is effectively eliminated in one relaxation 
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sweep and the convergence rate becomes that of the Pois- 
son equation for the pressure. It is possible to get ideal 
multigrid convergence rates for the system because each 
component of the error is treated appropriately. 

In the application of line Gauss-Seidel relaxation the 
unknowns along a given grid line are determined simulta- 
neously. A scalar tridiagonal matrix inversion is required 
for the pressure equation. Due to the second-order dis- 
cretization of the advection terms, a scalar pentadiagonal 
matrix inversion is necessary for the momentum equations. 

A straightforward FAS multigrid iteration is applied to 
the system of equations. Let La-i be the coarse grid op- 
erator, /*_! be the hne-to-coarse grid restriction operator, 
and /£ _1 be the coarse-to-hne grid prolongation operator. 
If qfc is the current solution on grid k, the residual on this 
grid is l'k = f k — La4a- This leads to the coarse-grid equa- 
tion 

La— iQa — 1 = fA — 1 = Ik - 1 r a + La— i {^k— i4a^ • (15) 

After solving the coarse-grid equation for Qa-i, the fine- 
grid solution is corrected by 

qA <- 4 a + /* _1 (qA-i - /a_i4a) • (16) 

Equation (15) is solved by applying the same relaxation 
procedure that is used to solve the fine-grid equation. 
Multigrid is applied recursively to the coarse-grid equa- 
tion. On the coarsest grid, many relaxation sweeps are 
performed to insure that the equation is solved completely. 
A conventional FF-cycle is used. 

Results 

The numerical scheme described in the previous sections 
has been applied to two incompressible, viscous flow prob- 
lems. As an initial evaluation of the scheme, high Reynolds 
number flow past a flat plate is considered. For the sec- 
ond problem flow past a symmetric Karman-Trefftz airfoil 
is solved on two different mesh topologies. 

The computational domain for the flat plate simulation 
is displayed in Fig. 2. The inflow boundary (ar = 0), outflow 
boundary (a~ = 3), and the upper boundary (y = 1) are lo- 
cated one plate length away from the plate. At the inflow 
boundary, the free-stream conditions (uoo = 1, = 0, 

Poo = —0.5) are specified. For the upper and downstream 
boundaries the pressure is set to p„o and the velocity com- 
ponents are obtained by solving the momentum equations. 
On the lower boundary symmetry conditions are applied 
upstream and downstream of the plate. The no-slip condi- 
tions are imposed on the plate (x = 1 to x = 2). A wake 
flow develops downstream of the plate that diffuses slowly 
in the present case of laminar flow. 

The finest grid used for the flat plate calculations con- 
sisted of 192 x 96 cells. In order to resolve the boundary 
layer on the plate the grid was clustered at y = 0 and 
stretched geometrically to y = 1 (see Fig. 2). For this grid 
the minimum spacing in the y direction is 0.002, and the 
stretching factor is 1.03. A series of nested coarse grids was 
obtained by coarsening the fine grids by a factor of two in 
each coordinate direction. In all cases shown below, the 
coarsest grid was 12 x 6 cells. 

Line Gauss-Seidel relaxation was used for the computa- 
tions. The complete relaxation process involves one sweep 



Fig. 2 Domain of flat plate flow with 48 x 24 cell 
grid. 



Fig. 3 Velocity profiles at midplate location for 
laminar flow over flat plate (Re = 10,000). 


with vertical line solves and one sweep with horizontal line 
solves, along with some additional work near the plate. 
The relatively small amount of additional work is needed 
due to the coupling of the flow equations near the plate. 
Moreover, the points for j < 3 are relaxed with about three 
sweeps of horizontal solves. For the vertical line solves 
relaxation begins at the inflow boundary (x = 0) and pro- 
ceeds to the outflow boundary (x = 3). In the case of the 
horizontal line solves relaxation starts at the outer bound- 
ary (y = 1) and continues to the inner boundary (y = 0). 
A If ( 1 . 1) multigrid cycle was used; that is, one complete 
relaxation process was performed on each grid before re- 
stricting to the coarse grid, and one complete relaxation 
process was performed after the coarse-grid correction was 
added to the fine-grid solution. 

The Reynolds number of the flow past the flat plate 
is 10, 000. Figure 3 shows the variation of the velocity 
component u nondimensionalized by the boundary-layer 
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17x17 33 x33 65x65 129 x129 257x257 513 x513 



Fig. 4 Convergence behavior for laminar flow over 
flat plate (Re = 10,000). Open symbols for pressure 
residuals; closed symbols for ^momentum residu- 
als. 



U/U e 


Fig. 5 Variation with multigrid cycles of velocity 
profile at midplate location for laminar flow over 
flat plate (Re = 10,000, 193 x 97 grid). 

edge velocity (« e ) with the scaled normal coordinate r]. 
The coordinate rj = 5y/S, where 5 is the thickness of the 
boundary layer. The computed velocity profiles are at the 
midplate location. Starting with the 96x48 cell grid, which 
has ten points in the boundary layer, there is excellent 
agreement with the classical Blasius solution. In Fig. 4 
the convergence behavior of the scheme is shown. The L 2 
norm of the residuals for the pressure and ^-momentum 
equations is given for each IT ( l . 1) cycle. For most of the 
grids the residuals have been reduced more than 6 orders of 
magnitude in 10 cycles. The convergence rate on the finest 



X 

Fig. 6 Near field of 128 x 64 cell O-type grid for 
computing laminar flow over Karman-Trefftz airfoil 
(.Re = 200). 

grid is 0.17 per cycle. As indicated in Fig. 5, the solution 
on the finest grid can nearly be obtained in a FMG process, 
involving only one multigrid cycle on each grid considered 
in Fig. 4. 

Solutions were obtained for incompressible, viscous flow 
(Re = 200 based on chord) around a nonlifting Karman- 
Trefftz airfoil. Previously, a form of the present numerical 
scheme was applied to inviscid airfoil flows, and an O-type 
mesh topology was used. So, initially in the present work, 
the same mesh topology was considered for viscous flows. 

A Karman-Trefftz airfoil was generated from a cylinder 
by a conformal mapping. 8 A trailing-edge angle of 10° is 
used, resulting in an airfoil of approximately 15% thick- 
ness. The airfoil flow is solved on a finite domain. At 
inflow points along the outer boundary the total pressure 
and flow inclination angle are specified. For outflow points 
the pressure is specified. The specified quantities are de- 
termined from the complex potential function for inviscid 
flow past the airfoil. 

A fine grid for the airfoil calculation (512 x 256 cells) 
was constructed by generating an O-grid in the circle plane 
and mapping it to the airfoil plane. The near field of the 
128 x 64 grid is displayed in Fig. 6. The outer boundary 
of the domain is roughly 3 chords from the airfoil. The 
coarsest grid in the grid sequence used for the multigrid 
solver contains 16 x 8 cells. On each grid a relaxation 
sweep started at the stagnation streamline, proceeded over 
the upper half of the domain, and then over the lower half 
of the domain. The relaxation sweep with radial solves was 
followed by sweeps with azimuthal solves near the airfoil, 
as in the flat plate case. 

Convergence behavior in the FMG process for the air- 
foil calculations with the O-type mesh is shown in Fig. 7. 
On the finest grid the asymptotic convergence rate for 
the x-momentum and pressure equations is about 0.14 per 
cycle. There is somewhat slower convergence on the coarser 
meshes, with about five orders of magnitude reduction of 
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Fig. 7 Convergence behavior for laminar flow over 
nonlifting Karman-Trefftz airfoil (Re = 200, O-type 
mesh topology). Open symbols for pressure resid- 
uals; closed symbols for ^-momentum residuals. 



X 

Fig. 8 Near field of 128 x 64 cell C-type grid for 
computing laminar flow over Karman-Trefftz airfoil 
(Re = 200). 

the residuals on the 64 x 32 cell grid. Such behavior was 
also observed by Roberts and Swanson 10 for inviscid flow. 

To facilitate the resolution of the near and far wake 
regions, especially in the case of laminar flow where the 
mixing rate due to the physical diffusion is slow, a C-type 
mesh was generated for the airfoil calculations. A hyper- 
bolic grid generator was used to create the mesh. The 
near field part of the grid is depicted in Fig. 8. In Fig. 9 
the convergence histories for this mesh topology are pre- 
sented. The average rate of reduction of the residual on 
the 512 x 256 grid is 0.19 for the pressure equation and 
0.22 for the ^-momentum equation. The slower rate for 
the ^-momentum equation is surprising since the rate of 
convergence of the scheme should be dictated by the ellip- 



Fig. 9 Convergence behavior for laminar flow over 
nonlifting Karman-Trefftz airfoil (Re = 200, C-type 
mesh topology). Open symbols for pressure resid- 
uals; closed symbols for i-momentum residuals. 

tic factor in this factorizable scheme. In addition, we do 
not observe this behavior with O-type mesh computation. 

In order to estimate the work required by the present 
and R-K schemes, we consider an operation count based on 
residual evaluations and the the number of cycles to achieve 
a specific reduction in residuals. For a W (m,n) multigrid 
cycle, there are m-fn residual evaluations on the fine grid. 
The restriction operator requires a residual calculation on 
the fine grid and one on the next coarser grid, with work 
corresponding to 1/4 that on the fine grid. In this assess- 
ment we neglect the cost of interpolating the residuals and 
solutions between the grid levels. For a W(m, n) cycle, we 
have that 9 

WU 1.11 

— — « m + n + 1+ -l l | - + - + •••) 

W- cycle v 4 2 4 ' 

= 2 (m + li + j), 

where WU means work unit. The WU per cycle for the cur- 
rent scheme with a W(2,l) cycle with only radial relaxation 
sweeps (neglecting additional work at surface boundary 
with azimuthal solves) is 8.5. With the five-stage R-K 
scheme and a W(1,0) cycle, the WU per cycle is 12.5. 

With the present scheme an additional residual evalu- 
ation is incurred due to the Gauss-Seidel updating of the 
solution. If we include that, the ratio of residual evalu- 
ations for the present scheme to that of the R-K scheme 
would be increased to 10.5 for the W(2,l) cycle used on 
the C-type grid and 8.5 for the W(l,l) cycle used for 
the O-type grid. Since we fully expect to have about the 
same number of residual evaluations for the different mesh 
topologies, the estimate of WU per cycle required by the 
present scheme is not altered. 

For the two schemes there is a difference in computa- 
tional effort resulting from scalar tridiagonal and pentadi- 
agonal inversions. As indicated previously each relaxation 
sweep of the present scheme requires one tridiagonal solve 
and two pentadiagonal solves (corresponding to the three 
flow equations) for the grid lines not in the direction of the 
relaxation. Since the operation count for a pentadiagonal 
inversion is about two times that for a tridiagonal inversion, 
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the work required on the fine grid with m-\-n = 3 is roughly 
proportional to that for 15 tridiagonal solves. For the R- 
K scheme there is a scalar tridiagonal inversion for each 
flow equation, each coordinate direction, and each of the 
five stages due to the residual smoothing procedure used 
to extend stability. So, in the particular case of three flow 
equations, the work due to this part of the scheme is pro- 
portional to 30 tridiagonal solves. Thus, the R-K scheme 
considered here 12 requires about two times as much work 
as the present scheme due to inversions. It should also 
be pointed out that the R-K scheme used here is a (5,3) 
scheme, which means that the complete residual, including 
viscous and numerical dissipation terms, is evaluated only 
on three stages. 

In order to roughly estimate the total speedup (i.e. , 
product of WU per cycle and number of cycles for desired 
level of convergence) of the present scheme relative to the 
R-K scheme, we not only consider residual reduction for 
the r-momentum equation but also the convergence of the 
drag coefficient (Cd). For a residual reduction of four or- 
ders of magnitude, the present scheme is at least an order 
magnitude faster than the R-K scheme for the laminar air- 
foil flow. This speedup is also evident in Figs. 10 and 11, 
which show the variation of Cd with multigrid cycles for the 
two schemes. The speedup of the present scheme on the 
finer meshes (mesh density > 256 x 128 cells) can exceed 
a factor of 25. 

Pressure and velocity contours for the solution on the 
512 x 256 C-type grid are depicted in Figs. 12 and 13. The 
thickness of the attached boundary layer is about 0.25 air- 
foil chords, and there are about 70 points in the boundary 
layer at the midchord location. Slow diffusion of the airfoil 
wake is evident. In Figs. 14 and 15 the pressure and skin 
friction computed with the present scheme on various grid 
densities are compared with the fine grid results obtained 
with the R-K scheme (where the free-stream Mach num- 
ber was set to 0.1). There is generally very good agreement 
starting with the 128 x 64 grid. 

Concluding Remarks 

A multigrid algorithm with essentially O(N) behavior 
has been developed for the steady Navier-Stokes equations. 
It has the virtue of simplicity; conventional finite-difference 
or finite-volume discretizations of the governing equations 
may be used, allowing flexibility in the choice of the under- 
lying numerical method. Appropriate discretization and 
efficient treatment of the pressure boundary condition has 
been demonstrated. Solutions have been obtained for lami- 
nar flow past a flat plate with essentially texbook multigrid 
efficiency. Furthermore, the multigrid method has been 
used to solve viscous flow with curvature effects. Fast con- 
vergence was obtained for this case also. With the present 
scheme more than an order of magnitude reduction in com- 
putational effort has been achieved when compared with a 
R-K based multigrid scheme. 

Further studies of the present scheme are required to de- 
termine the effects of mesh aspect ratio and grid stretching 
on convergence rate, and thus, on the discrete factorizabil- 
ity of the algorithm. While significant gains in computa- 
tional efficiency have been achieved for laminar flow, much 
greater gains (i.e., two orders of magnitude) are anticipated 


with the extension of factoriazable schemes to allow com- 
putation of turbulent flows. 
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Fig. 10 Convergence history of total drag for F! S' , 12 Pressure contours for laminar flow 

, . a - - rr m f nonlifting Karman-Trefftz airfoil (Re = 200). 

laminar flow over nonlitting Karman- iretttz airtoil ° v ' 

(Re = 200). 
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Fig. 11 Convergence history with R-K scheme of 
total drag for laminar flow over nonlifting Karman- 
Trefftz airfoil (Re = 200). 


Fig. 13 Velocity contours for laminar flow 
nonlifting Karman-Trefftz airfoil (Re = 200). 
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Fig. 14 Surface pressure distribution for laminar 
flow over nonlifting Karman-Trefftz airfoil (Re = 

200 ) . 


Karman-Trefftz Airfoil, Re c = 200, a = 0° 
Surface Skin Friction 
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Fig. 15 Surface skin-friction distribution for lam- 
inar flow over nonlifting Karman-Trefftz airfoil 
(Re = 200). 
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